using System;
using System.Diagnostics;
using System.Collections.Generic;

namespace PortfolioAllocator.Optimization
{
    /// <summary>Find solution of a (convex) Quadratic Programming problem</summary>
    class QPSolver
    {
        private static readonly double EPSILON = 1e-16;

        /// <summary>Implements the algorithm of Goldfarb and Idnani for the solution of a (convex) Quadratic Programming problem.
        /// The problem is in the form:
        /// min 0.5 * x G x + g0 x
        /// s.t.
        /// CE^T x + ce0 = 0
        /// CI^T x + ci0 >= 0
        /// The matrix and vectors dimensions are as follows:
        /// G: n * n
        /// g0: n				
        /// CE: n * p
        /// ce0: p
        ///	CI: n * m
        /// ci0: m
        /// x: n
        /// </summary>
        /// <remark>
        /// Notes:
        /// 1. pay attention in setting up the vectors ce0 and ci0. 
        /// If the constraints of your problem are specified in the form 
        /// A^T x = b and C^T x >= d, then you should set ce0 = -b and ci0 = -d.
        /// 2. The matrix G is modified within the function since it is used to compute
        /// the G = L^T L cholesky factorization for further computations inside the function. 
        /// If you need the original matrix G you should make a copy of it and pass the copy
        /// to the function.
        /// </remark>        
        /// <param name="G">n*n Quadratic part of the objective function</param>
        /// <param name="g0">n dimensional linear part of the objective function</param>
        /// <param name="CE">n*p matrix from Equality constraint</param>
        /// <param name="ce0">p dimensional constant part of equality constraint</param>
        /// <param name="CI">n*m matrix form in-equality constraint</param>
        /// <param name="ci0">m dimensional constant part of inequality constraint</param>
        /// <param name="x">n dimensional result</param>
        /// <returns>
        /// The function will return the cost of the solution written in the x vector or
        /// std::numeric_limits::infinity() if the problem is infeasible. In the latter case
        /// the value of the x vector is not correct.
        /// </returns>
        public double solve(double[,] G, double[] g0,
                double[,] CE, double[] ce0,
                double[,] CI, double[] ci0,
                out double[] x)
        {

            int n = G.GetLength(0);
            int p = CE.GetLength(1);
            int m = CI.GetLength(1);
            Debug.Assert(G.GetLength(0) == G.GetLength(1), "G is not square matrix", "G(n,n) = G(" + G.GetLength(0) + "," + G.GetLength(1) + ")");
            Debug.Assert(g0.Length == n, "g.Length is not equal to G dimension", "g0.Length=" + g0.Length + ", expected n=" + n);
            Debug.Assert(CE.GetLength(0) == n, "CE number of rows is not equal to G dimension", "CE(n,p)=CE(" + CE.GetLength(0) + "," + CE.GetLength(1) + "). Expected CE(n=" + n + ",p=" + p + ")");
            Debug.Assert(ce0.Length == p, "ce.Length is not equal to CE number of columns", "ce(p)=ce(" + ce0.Length + "). Expected ce(p=" + p + ")");
            Debug.Assert(CI.GetLength(0) == n, "CI number of rows is not equal to G dimension", "CI(n,m)=CI(" + CI.GetLength(0) + "," + CI.GetLength(1) + "). Expected CI(n=" + n + ",m=" + m + ")");
            Debug.Assert(ci0.Length == m, "ci.Length is not equal to CE number of columns", "ci(m)=ci(" + ci0.Length + "). Expected ci(m=" + m + ")");

            int ip, me, mi;
            x = new double[n];  // result
            double[,] R = new double[n, n];

            double[] s = new double[m + p];
            double[] z = new double[n];
            double[] r = new double[m + p];
            double[] d = new double[n];

            double[,] J = new double[n, n];

            double[] np = new double[n];
            double[] u = new double[m + p];
            double[] x_old = new double[n];
            double[] u_old = new double[m + p];
            double psi, sum, ss, R_norm, c1, c2;
            const double inf = double.PositiveInfinity;
            double t, t1, t2; /* t is the step lenght, which is the minimum of the partial step length t1 
					  * and the full step length t2 */
            int[] A = new int[m + p];
            int[] A_old = new int[m + p];
            int[] iai = new int[m + p];
            int q;
            int iq;
            bool[] iaexcl = new bool[m + p];


            me = p; /* number of equality constraints */
            mi = m; /* number of inequality constraints */
            q = 0;  /* size of the active set A (containing the indices of the active constraints) */

            preprocessProblem(G, J, R, d, z, out c1, out c2, out R_norm);
            double f_value = findUnconstrainedSolution(G, x, g0);
            addEqualityConstraintsToWorkingSet(n, me, mi, CE, J, R, ce0, x, d, z, u, iai, A, np, r, out t2, ref f_value, ref R_norm, out iq);

            // The main iteration loop (l1)
            const int maxIter = 1000;
            for (int iter = 0; iter < maxIter; ++iter)
            {
                /* step 1: choose a violated constraint */
                for (int i = me; i < iq; i++)
                {
                    ip = A[i];
                    iai[ip] = -1;
                }

                /* compute s[x] = ci^T * x + ci0 for all elements of K \ A */
                ss = 0.0;
                psi = 0.0; /* this value will contain the sum of all infeasibilities */
                ip = 0; /* ip will be the index of the chosen violated constraint */
                for (int i = 0; i < mi; i++)
                {
                    iaexcl[i] = true;
                    sum = 0.0;
                    for (int j = 0; j < n; j++)
                        sum += CI[j, i] * x[j];
                    sum += ci0[i];
                    s[i] = sum;
                    psi += Math.Min(0.0, sum);
                }

                if (Math.Abs(psi) <= mi * double.Epsilon * c1 * c2 * 100.0)
                // numerically there are not infeasibilities anymore 
                {
                    q = iq;
                    return f_value;
                }

                /* save old values for u and A */
                for (int i = 0; i < iq; i++)
                {
                    u_old[i] = u[i];
                    A_old[i] = A[i];
                }
                /* and for x */
                for (int i = 0; i < n; i++)
                    x_old[i] = x[i];
                /* Step 2: check for feasibility and determine a new S-pair */
                while (true) // (l2) loop
                {
                    for (int i = 0; i < mi; i++)
                    {
                        if (s[i] < ss && iai[i] != -1 && iaexcl[i])
                        {
                            ss = s[i];
                            ip = i;
                        }
                    }
                    if (ss >= 0.0)
                    {
                        q = iq;
                        return f_value;
                    }

                    /* set np = n[ip] */
                    for (int i = 0; i < n; i++)
                        np[i] = CI[i, ip];
                    /* set u = [u 0]^T */
                    u[iq] = 0.0;
                    /* add ip to the active set A */
                    A[iq] = ip;

                    /* Step 2a: determine step direction */
                    /* compute z = H np: the step direction in the primal space (through J, see the paper) */
                    while (true) // (l2a) loop
                    {
                        compute_d(d, J, np, n);
                        update_z(z, J, d, n, iq);
                        /* compute N* np (if q > 0): the negative of the step direction in the dual space */
                        update_r(R, r, d, n, iq);
                        /* Step 2b: compute step length */
                        int l = 0;
                        /* Compute t1: partial step length (maximum step in dual space without violating dual feasibility */
                        t1 = inf; /* +inf */
                        /* find the index l s.t. it reaches the minimum of u+[x] / r */
                        for (int k = me; k < iq; k++)
                        {
                            if (r[k] > 0.0)
                            {
                                if (u[k] / r[k] < t1)
                                {
                                    t1 = u[k] / r[k];
                                    l = A[k];
                                }
                            }
                        }
                        /* Compute t2: full step length (minimum step in primal space such that the constraint ip becomes feasible */
                        if (Math.Abs(scalar_product(z, z, n)) > double.Epsilon) // i.e. z != 0
                            t2 = -s[ip] / scalar_product(z, np, n);
                        else
                            t2 = inf; /* +inf */

                        /* the step is chosen as the minimum of t1 and t2 */
                        t = Math.Min(t1, t2);
                        /* Step 2c: determine new S-pair and take step: */
                        /* case (i): no step in primal or dual space */
                        if (t >= inf)
                        {
                            /* QPP is infeasible */
                            // unbounded to raise
                            q = iq;
                            // return inf;
                            throw new Exception("Infeasible problem");
                        }

                        /* case (ii): step in dual space */
                        if (t2 >= inf)
                        {
                            /* set u = u +  t * [-r 1] and drop constraint l from the active set A */
                            for (int k = 0; k < iq; k++)
                                u[k] -= t * r[k];
                            u[iq] += t;
                            iai[l] = l;
                            delete_constraint(R, J, A, u, n, p, ref iq, l);
                        } // end of 'if(t2 >= inf)'
                        else // t2 < inf
                        {
                            /* case (iii): step in primal and dual space */
                            /* set x = x + t * z */
                            for (int k = 0; k < n; k++)
                                x[k] += t * z[k];
                            /* update the solution value */
                            f_value += t * scalar_product(z, np, n) * (0.5 * t + u[iq]);
                            /* u = u + t * [-r 1] */
                            for (int k = 0; k < iq; k++)
                                u[k] -= t * r[k];
                            u[iq] += t;

                            if (t == t2)
                                break; // break from l2a loop

                            /* a patial step has taken */
                            /* drop constraint l */
                            iai[l] = l;
                            delete_constraint(R, J, A, u, n, p, ref iq, l);
                            /* update s[ip] = CI * x + ci0 */
                            sum = 0.0;
                            for (int k = 0; k < n; k++)
                                sum += CI[k, ip] * x[k];
                            s[ip] = sum + ci0[ip];
                        }// end of else t2 < inf
                    } // end of l2a loop (goto l2a)	

                    //here (t == t2 < inf)
                    /* full step has taken */
                    /* add constraint ip to the active set*/
                    if (add_constraint(R, J, d, n, ref iq, ref R_norm))
                    {
                        iai[ip] = -1;
                        break; // go to loop l1
                    }

                    // we failed to add constraint -> continue with l2
                    iaexcl[ip] = false;
                    delete_constraint(R, J, A, u, n, p, ref iq, ip);
                    for (int i = 0; i < m; i++)
                        iai[i] = i;
                    for (int i = 0; i < iq; i++)
                    {
                        A[i] = A_old[i];
                        iai[A[i]] = -1;
                        u[i] = u_old[i];
                    }
                    for (int i = 0; i < n; i++)
                        x[i] = x_old[i];
                } // endl of l2 loop (goto l2)
            }// end of iter loop (goto l1)
            throw new Exception("Max number of internal iteration reached without result");


        }

        #region Private Methods
        /// <summary>Preprocessing phase</summary>
        /// <param name="G">defines the quadratic term in QP. out -> cholesky decomposed G</param>
        /// <param name="J">inverse of factorized G^-1</param>
        /// <param name="c1">c1 outputs trace of the original matrix G</param>
        private void preprocessProblem(double[,] G, double[,] J, double[,] R,
                               double[] d, double[] z,
                               out double c1, out double c2, out double R_norm)
        {
            int n = G.GetLength(0);
            /* compute the trace of the original matrix G */
            c1 = 0.0;
            for (int i = 0; i < n; i++)
                c1 += G[i, i];

            /* decompose the matrix G in the form L^T L */
            cholesky_decomposition(G, n);

            /* initialize the matrix R */
            for (int i = 0; i < n; i++)
            {
                d[i] = 0.0;
                for (int j = 0; j < n; j++)
                    R[i, j] = 0.0;
            }
            R_norm = 1.0; /* this variable will hold the norm of the matrix R */

            /* compute the inverse of the factorized matrix G^-1, this is the initial value for H */
            c2 = 0.0;
            for (int i = 0; i < n; i++)
            {
                d[i] = 1.0;
                forward_elimination(G, z, d, n);
                for (int j = 0; j < n; j++)
                    J[i, j] = z[j];
                c2 += z[i];
                d[i] = 0.0;
            }
            /* c1 * c2 is an estimate for cond(G) */
        }

        /// <summary>
        /// Find the unconstrained minimizer of the quadratic form 0.5 * x G x + g0 x  
        /// this is a feasible point in the dual space x = G^-1 * g0
        /// </summary>
        /// <returns>unconstrained solution value</returns>
        private double findUnconstrainedSolution(double[,] G, double[] x, double[] g0)
        {
            int n = G.GetLength(0);
            cholesky_solve(G, x, g0, n);
            for (int i = 0; i < n; i++)
                x[i] = -x[i];
            /* and compute the current solution value */
            double f_value = 0.5 * scalar_product(g0, x, n);
            return f_value;
        }

        /// <summary>Add equality constraints to the working set A </summary>
        private void addEqualityConstraintsToWorkingSet(int n, int me, int mi,
                                           double[,] CE, double[,] J, double[,] R,
                                           double[] ce0,
                                           double[] x, double[] d, double[] z, double[] u,
                                           int[] iai, int[] A, double[] np, double[] r,
                                           out double t2, ref double f_value, ref double R_norm, out int iq)
        {
            iq = 0;
            t2 = 0.0;
            for (int i = 0; i < me; i++)
            {
                for (int j = 0; j < n; j++)
                    np[j] = CE[j, i];
                compute_d(d, J, np, n);
                update_z(z, J, d, n, iq);
                update_r(R, r, d, n, iq);

                /* compute full step length t2: i.e., the minimum step in primal space s.t. the contraint 
                becomes feasible */
                t2 = 0.0;
                if (Math.Abs(scalar_product(z, z, n)) > double.Epsilon) // i.e. z != 0
                    t2 = (-scalar_product(np, x, n) - ce0[i]) / scalar_product(z, np, n);

                /* set x = x + t2 * z */
                for (int k = 0; k < n; k++)
                    x[k] += t2 * z[k];

                /* set u = u+ */
                u[iq] = t2;
                for (int k = 0; k < iq; k++)
                    u[k] -= t2 * r[k];

                /* compute the new solution value */
                f_value += 0.5 * (t2 * t2) * scalar_product(z, np, n);
                A[i] = -i - 1;

                if (!add_constraint(R, J, d, n, ref iq, ref R_norm))
                    throw new Exception("Equality constraints are linearly dependent");
            }

            /* set iai = K \ A */
            for (int i = 0; i < mi; i++)
                iai[i] = i;
        }


        private void compute_d(double[] d, double[,] J, double[] np, int n)
        {
            int i, j;
            double sum;

            /* compute d = H^T * np */
            for (i = 0; i < n; i++)
            {
                sum = 0.0;
                for (j = 0; j < n; j++)
                    sum += J[j, i] * np[j];
                d[i] = sum;
            }
        }

        private void update_z(double[] z, double[,] J, double[] d, int n, int iq)
        {
            int i, j;

            /* setting of z = H * d */
            for (i = 0; i < n; i++)
            {
                z[i] = 0.0;
                for (j = iq; j < n; j++)
                    z[i] += J[i, j] * d[j];
            }
        }

        private void update_r(double[,] R, double[] r, double[] d, int n, int iq)
        {
            int i, j;
            double sum;

            /* setting of r = R^-1 d */
            for (i = iq - 1; i >= 0; i--)
            {
                sum = 0.0;
                for (j = i + 1; j < iq; j++)
                    sum += R[i, j] * r[j];
                r[i] = (d[i] - sum) / R[i, i];
            }
        }

        private bool add_constraint(double[,] R, double[,] J, double[] d, int n, ref int iq, ref double R_norm)
        {
            int i, j, k;
            double cc, ss, h, t1, t2, xny;

            /* we have to find the Givens rotation which will reduce the element
            d[j] to zero.
            if it is already zero we don't have to do anything, except of
            decreasing j */
            for (j = n - 1; j >= iq + 1; j--)
            {
                /* The Givens rotation is done with the matrix (cc cs, cs -cc).
                If cc is one, then element (j) of d is zero compared with element
                (j - 1). Hence we don't have to do anything. 
                If cc is zero, then we just have to switch column (j) and column (j - 1) 
                of J. Since we only switch columns in J, we have to be careful how we
                update d depending on the sign of gs.
                Otherwise we have to apply the Givens rotation to these columns.
                The i - 1 element of d has to be updated to h. */
                cc = d[j - 1];
                ss = d[j];
                h = distance(cc, ss);
                if (h == 0.0)
                    continue;
                d[j] = 0.0;
                ss = ss / h;
                cc = cc / h;
                if (cc < 0.0)
                {
                    cc = -cc;
                    ss = -ss;
                    d[j - 1] = -h;
                }
                else
                    d[j - 1] = h;
                xny = ss / (1.0 + cc);
                for (k = 0; k < n; k++)
                {
                    t1 = J[k, j - 1];
                    t2 = J[k, j];
                    J[k, j - 1] = t1 * cc + t2 * ss;
                    J[k, j] = xny * (t1 + J[k, j - 1]) - t2;
                }
            }
            /* update the number of constraints added*/
            iq++;
            /* To update R we have to put the iq components of the d vector
            into column iq - 1 of R
            */
            for (i = 0; i < iq; i++)
                R[i, iq - 1] = d[i];
            if (Math.Abs(d[iq - 1]) <= EPSILON * R_norm)  //  if (Math.Abs(d[iq - 1]) <= Double.Epsilon * R_norm)
                // problem degenerate
                return false;
            R_norm = Math.Max(R_norm, Math.Abs(d[iq - 1]));
            return true;
        }

        private void delete_constraint(double[,] R, double[,] J, int[] A, double[] u, int n, int p, ref int iq, int l)
        {

            int i, j, k, qq = 0;
            double cc, ss, h, xny, t1, t2;

            /* Find the index qq for active constraint l to be removed */
            for (i = p; i < iq; i++)
                if (A[i] == l)
                {
                    qq = i;
                    break;
                }

            /* remove the constraint from the active set and the duals */
            for (i = qq; i < iq - 1; i++)
            {
                A[i] = A[i + 1];
                u[i] = u[i + 1];
                for (j = 0; j < n; j++)
                    R[j, i] = R[j, i + 1];
            }

            A[iq - 1] = A[iq];
            u[iq - 1] = u[iq];
            A[iq] = 0;
            u[iq] = 0.0;
            for (j = 0; j < iq; j++)
                R[j, iq - 1] = 0.0;
            /* constraint has been fully removed */
            iq--;

            if (iq == 0)
                return;

            for (j = qq; j < iq; j++)
            {
                cc = R[j, j];
                ss = R[j + 1, j];
                h = distance(cc, ss);
                if (h == 0.0)
                    continue;
                cc = cc / h;
                ss = ss / h;
                R[j + 1, j] = 0.0;
                if (cc < 0.0)
                {
                    R[j, j] = -h;
                    cc = -cc;
                    ss = -ss;
                }
                else
                    R[j, j] = h;

                xny = ss / (1.0 + cc);
                for (k = j + 1; k < iq; k++)
                {
                    t1 = R[j, k];
                    t2 = R[j + 1, k];
                    R[j, k] = t1 * cc + t2 * ss;
                    R[j + 1, k] = xny * (t1 + R[j, k]) - t2;
                }
                for (k = 0; k < n; k++)
                {
                    t1 = J[k, j];
                    t2 = J[k, j + 1];
                    J[k, j] = t1 * cc + t2 * ss;
                    J[k, j + 1] = xny * (J[k, j] + t1) - t2;
                }
            }
        }

        private double distance(double a, double b)
        {
            double a1, b1, t;
            a1 = Math.Abs(a);
            b1 = Math.Abs(b);
            if (a1 > b1)
            {
                t = (b1 / a1);
                return a1 * Math.Sqrt(1.0 + t * t);
            }
            else
                if (b1 > a1)
                {
                    t = (a1 / b1);
                    return b1 * Math.Sqrt(1.0 + t * t);
                }
            return a1 * Math.Sqrt(2.0);
        }

        private double scalar_product(double[] x, double[] y, int n)
        {
            int i;
            double sum;

            sum = 0.0;
            for (i = 0; i < n; i++)
                sum += x[i] * y[i];
            return sum;
        }

        private void cholesky_decomposition(double[,] A, int n)
        {
            double[] d = new double[n];
            int i, j, k;
            double sum;

            for (i = 0; i < n; i++)
            {
                for (j = i; j < n; j++)
                {
                    sum = A[i, j];
                    for (k = i - 1; k >= 0; k--)
                        sum -= A[i, k] * A[j, k];
                    if (i == j)
                    {
                        if (sum <= 0.0)
                        {
                            // raise error
                            throw new Exception("The matrix passed to the Cholesky A = L L^T decomposition is not positive definite");
                        }
                        A[i, i] = Math.Sqrt(sum);
                    }
                    else
                        A[j, i] = sum / A[i, i];
                }
                for (k = i + 1; k < n; k++)
                    A[i, k] = A[k, i];
            }
        }

        private void cholesky_solve(double[,] L, double[] x, double[] b, int n)
        {
            double[] y = new double[n];

            /* Solve L * y = b */
            forward_elimination(L, y, b, n);
            /* Solve L^T * x = y */
            backward_elimination(L, x, y, n);
        }

        private void forward_elimination(double[,] L, double[] y, double[] b, int n)
        {
            int i, j;

            y[0] = b[0] / L[0, 0];
            for (i = 1; i < n; i++)
            {
                y[i] = b[i];
                for (j = 0; j < i; j++)
                    y[i] -= L[i, j] * y[j];
                y[i] = y[i] / L[i, i];
            }
        }

        private void backward_elimination(double[,] U, double[] x, double[] y, int n)
        {
            int i, j;

            x[n - 1] = y[n - 1] / U[n - 1, n - 1];
            for (i = n - 2; i >= 0; i--)
            {
                x[i] = y[i];
                for (j = i + 1; j < n; j++)
                    x[i] -= U[i, j] * x[j];
                x[i] = x[i] / U[i, i];
            }
        }
        #endregion
    }

    /// <summary>Definition of position constraint:  max >= sum x_i >= min</summary>
    struct PositionConstraint
    {
        private double minM;
        private double maxM;
        private List<int> ixM;  // list of indexes of positions in the constraint

        public PositionConstraint(double min, double max, List<int> indexes)
        {
            this.minM = min;
            this.maxM = max;
            this.ixM = indexes;
        }
        public PositionConstraint(double min, double max, int[] indexes)
            : this(min, max, new List<int>(indexes))
        { }

        public double Min
        {
            get { return this.minM; }
        }
        public double Max
        {
            get { return this.maxM; }
        }
        public List<int> Indexes
        {
            get { return this.ixM; }
        }
    }

    class AllocationSolver
    {
        /// <summary>Find portfolio weights minimizing variance given Position Constraints
        ///     min:     1/2 x^T cov x      
        ///     subj to: xMax >= sum_i x_i >= xMin 
        ///              sum_i x_i = 1 
        /// Returns an array with the optimal portfolio weights.
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="constraints">List of position constraints</param>
        /// <returns>vector of optimal weights</returns>
        public double[] findMinVariancePortfWeights(double[,] cov, List<PositionConstraint> constraints)
        {
            int n = cov.GetLength(0);
            if (n != cov.GetLength(1))
                throw new Exception("Inconsistent input matrix");

            // 1. objective function f(x) = 1/2 xGx + 0x
            double[,] G = (double[,])cov.Clone();  // create a writeable copy
            double[] g = new double[n]; // obj function linear term = 0	

            // 2. Equality constraints sum_i x_i = 1    ( sum_i 1*x_i - 1 = 0)
            double[,] CE = new double[n, 1];
            for (int i = 0; i < n; ++i)
                CE[i, 0] = 1;
            double[] ce = { -1 };

            // 3. InEquality constraints xMin <= sum_indexes x_i <= xMax 
            int p = 2 * constraints.Count;
            double[,] CI = new double[n, p];
            double[] ci = new double[p];
            for (int i = 0; i < constraints.Count; ++i)
            {
                PositionConstraint c = constraints[i];
                
                ci[2 * i] = -c.Min;
                ci[2 * i + 1] = c.Max;
                foreach (int ix in c.Indexes)
                {
                    CI[ix, 2 * i] = +1;         // sum_indexes (+1 * x_i) - xMin >= 0
                    CI[ix, 2 * i + 1] = -1;     // sum_indexes (-1 * x_i) + Max >= 0 
                }
            }

            double[] x;
            QPSolver solver = new QPSolver();
            double f = solver.solve(G, g, CE, ce, CI, ci, out x);
            return x;
        }

        /// <summary>Find portfolio weights minimizing variance:
        ///     min:     1/2 x^T cov x      
        ///     subj to: xMax_i >= x_i >= xMin_i 
        ///              sum_i x_i = 1 
        /// Returns an array with the optimal portfolio weights.
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="xMin">vector of minimum weights</param>
        /// <param name="xMax">vector of maximum weights</param>
        /// <returns>vector of optimal weights</returns>
        public double[] findMinVariancePortfWeights(double[,] cov, double[] xMin, double[] xMax)
        {
            return this.findMinVariancePortfWeights(cov, null, 0, xMin, xMax);
        }

        /// <summary>Find portfolio weights minimizing variance given level of expected return:
        ///     min:     1/2 x^T cov x      
        ///     subj to: x^T r >= rMin
        ///              xMax_i >= x_i >= xMin_i 
        ///              sum_i x_i = 1 
        /// Returns an array with the optimal portfolio weights.
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="r">Vector of expected returns. If r==null then return constraint is ignored</param>
        /// <param name="rMin">minimum level of expected return</param>
        /// <param name="xMin">vector of minimum weights</param>
        /// <param name="xMax">vector of maximum weights</param>
        /// <returns>vector of optimal weights</returns>
        public double[] findMinVariancePortfWeights(double[,] cov, double[] r, double rMin, double[] xMin, double[] xMax)
        {
            bool useRetConstr = (r != null);
            int n = cov.GetLength(0);
            if (n != cov.GetLength(1) || n != xMin.Length || n != xMax.Length || (useRetConstr && n != r.Length))
                throw new Exception("Inconsistent input matrix / vectors dimensions");

            // 1. objective function f(x) = 1/2 xGx + 0x
            double[,] G = (double[,]) cov.Clone();  // create a writeable copy
            double[] g = new double[n]; // obj function linear term = 0	

            // 2. Equality constraints sum_i x_i = 1    ( sum_i 1*x_i - 1 = 0)
            double[,] CE = new double[n, 1];
            for (int i = 0; i < n; ++i)
                CE[i, 0] = 1;
            double[] ce = { -1 };

            // 3. InEquality constraints xMin <= x <= xMax and x r >= rMin (if r!=null)
            int p = (useRetConstr ? 2 * n + 1 : 2 * n);
            double[,] CI = new double[n, p];
            double[] ci = new double[p];
            for (int i = 0; i < n; ++i)
            {
                // x[i] - xMin[i] >= 0		(xMin[i] <= x[i])
                CI[i, 2 * i] = 1;
                ci[2 * i] = -xMin[i];
                // -x[i] + xMax[i] >= 0		(x[i] <= xMax[i])
                CI[i, 2 * i + 1] = -1;
                ci[2 * i + 1] = xMax[i];
                if (useRetConstr)
                    // sum_i r[i] x[i] >= rMin
                    CI[i, p - 1] = r[i];
            }
            if (useRetConstr)
                ci[p - 1] = -rMin;

            double[] x;
            QPSolver solver = new QPSolver();
            double f = solver.solve(G, g, CE, ce, CI, ci, out x);
            return x;
        }

        /// <summary>Find portfolio weights maximizing the expected log return (Kelly optimal portfolio):
        ///     max: x * r - 1/2 x^T cov x
        ///     subj to:    xMax_i >= x_i >= xMin_i
        ///                 sum_i x_i = 1
        /// Returns an array with the optimal portfolio weights.
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="r">Vector of expected returns</param>
        /// <param name="xMin">vector of minimum weights</param>
        /// <param name="xMax">vector of maximum weights</param>
        /// <returns>vector of optimal weights</returns>
        public double[] findMaxLogReturnPortfWeights(double[,] cov, double[] r, List<PositionConstraint> constraints)
        {
            int n = cov.GetLength(0);
            if (n != cov.GetLength(1) || n != r.Length )
                throw new Exception("Inconsistent input matrix / vectors dimensions");

            // 1. objective function f(x) = 1/2 xGx - rx
            double[,] G = (double[,])cov.Clone();  // create a writeable copy
            double[] g = new double[n]; // obj function linear term: -r*x	
            for (int i = 0; i < n; ++i)
                g[i] = -r[i];

            // 2. Equality constraints sum_i x_i = 1    ( sum_i 1*x_i - 1 = 0)
            double[,] CE = new double[n, 1];
            for (int i = 0; i < n; ++i)
                CE[i, 0] = 1;
            double[] ce = { -1 };

            // 3. InEquality constraints xMin <= sum_indexes x_i <= xMax 
            int p = 2 * constraints.Count;
            double[,] CI = new double[n, p];
            double[] ci = new double[p];
            for (int i = 0; i < constraints.Count; ++i)
            {
                PositionConstraint c = constraints[i];

                ci[2 * i] = -c.Min;
                ci[2 * i + 1] = c.Max;
                foreach (int ix in c.Indexes)
                {
                    CI[ix, 2 * i] = +1;         // sum_indexes (+1 * x_i) - xMin >= 0
                    CI[ix, 2 * i + 1] = -1;     // sum_indexes (-1 * x_i) + Max >= 0 
                }
            }

            double[] x;
            QPSolver solver = new QPSolver();
            double f = solver.solve(G, g, CE, ce, CI, ci, out x);
            return x;
        }

        /// <summary>Find portfolio weights maximizing the expected log return (Kelly optimal portfolio):
        ///     max: x * r - 1/2 x^T cov x
        ///     subj to:    xMax_i >= x_i >= xMin_i
        ///                 sum_i x_i = 1
        /// Returns an array with the optimal portfolio weights.
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="r">Vector of expected returns</param>
        /// <param name="xMin">vector of minimum weights</param>
        /// <param name="xMax">vector of maximum weights</param>
        /// <returns>vector of optimal weights</returns>
        public double[] findMaxLogReturnPortfWeights(double[,] cov, double[] r, double[] xMin, double[] xMax)
        {
            int n = cov.GetLength(0);
            if (n != cov.GetLength(1) || n != r.Length || n != xMin.Length || n != xMax.Length)
                throw new Exception("Inconsistent input matrix / vectors dimensions");

            // 1. objective function f(x) = 1/2 xGx - rx
            double[,] G = (double[,])cov.Clone();  // create a writeable copy
            double[] g = new double[n]; // obj function linear term: -r*x	
            for (int i = 0; i < n; ++i)
                g[i] = -r[i];

            // 2. Equality constraints sum_i x_i = 1    ( sum_i 1*x_i - 1 = 0)
            double[,] CE = new double[n, 1];
            for (int i = 0; i < n; ++i)
                CE[i, 0] = 1;
            double[] ce = { -1 };

            // 3. InEquality constraints xMin <= x <= xMax
            int p = 2 * n;
            double[,] CI = new double[n, p];
            double[] ci = new double[p];
            for (int i = 0; i < n; ++i)
            {
                // x[i] - xMin[i] >= 0		(xMin[i] <= x[i])
                CI[i, 2 * i] = 1;
                ci[2 * i] = -xMin[i];
                // -x[i] + xMax[i] >= 0		(x[i] <= xMax[i])
                CI[i, 2 * i + 1] = -1;
                ci[2 * i + 1] = xMax[i];
            }

            double[] x;
            QPSolver solver = new QPSolver();
            double f = solver.solve(G, g, CE, ce, CI, ci, out x);
            return x;
        }

        /// <summary>Find portfolio weights INCLUDING RISK FREE ASSET maximizing the expected log return (Kelly optimal portfolio):
        ///     max: (1-x)*rFree + x * r - 1/2 x^T cov x
        ///     subj to:    xMax_i >= x_i >= xMin_i
        /// Returns an array with the optimal portfolio weights. Risk free wights corresponds to 1 - sum(x).
        /// </summary>
        /// <param name="cov">Covariance matrix</param>
        /// <param name="r">Vector of expected returns</param>
        /// <param name="rFree">Risk free return</param>
        /// <param name="xMin">vector of minimum weights</param>
        /// <param name="xMax">vector of maximum weights</param>
        /// <returns>vector of optimal weights on risky assets. Risk-free weight = 1 - sum(x)</returns>
        public double[] findMaxLogReturnPortfWeights(double[,] cov, double[] r, double rFree, double[] xMin, double[] xMax)
        {
            int n = cov.GetLength(0);
            if (n != cov.GetLength(1) || n != r.Length || n != xMin.Length || n != xMax.Length)
                throw new Exception("Inconsistent input matrix / vectors dimensions");

            // 1. objective function f(x) = 1/2 xGx - rx - rFree*(1-x) = 1/2 xGx - (r-rFree)x - rFree
            double[,] G = (double[,])cov.Clone();  // create a writeable copy
            double[] g = new double[n]; // obj function linear term: -r*x	
            for (int i = 0; i < n; ++i)
                g[i] = -(r[i] - rFree);

            double[,] CE = new double[n, 0];
            double[] ce = new double[0];

            // 2. InEquality constraints xMin <= x <= xMax
            int p = 2 * n;
            double[,] CI = new double[n, p];
            double[] ci = new double[p];
            for (int i = 0; i < n; ++i)
            {
                // x[i] - xMin[i] >= 0		(xMin[i] <= x[i])
                CI[i, 2 * i] = 1;
                ci[2 * i] = -xMin[i];
                // -x[i] + xMax[i] >= 0		(x[i] <= xMax[i])
                CI[i, 2 * i + 1] = -1;
                ci[2 * i + 1] = xMax[i];
            }

            double[] x;
            QPSolver solver = new QPSolver();
            double f = solver.solve(G, g, CE, ce, CI, ci, out x);
            return x;
        }

        /// <summary>Calculate portfolio expected return from portfolio weights and instrument expected returns</summary>
        /// <param name="x">Portfolio instrument weights</param>
        /// <param name="r">Instrument expected returns</param>
        /// <returns>Portoflio expected return = sum_i x_i * r_i</returns>
        public double calcExpReturn(double[] x, double[] r)
        {
            int n = x.Length;
            if (n != r.Length)
                throw new System.IndexOutOfRangeException("Weights and expected returns have different sizes");
            double expRet = 0;
            for (int i = 0; i < n; ++i)
                expRet += x[i] * r[i];
            return expRet;
        }

        /// <summary>Calculate portfolio expected volatility from portfolio weights and covariance matrix</summary>
        /// <param name="x">Portfolio instrument weights</param>
        /// <param name="cov">Covariance matrix</param>
        /// <returns>Portoflio expected volatility = sum_ij x_i * cov_ij * x_j </returns>
        public double calcVolatility(double[] x, double[,] cov)
        {
            int n = x.Length;
            if (n != cov.GetLength(0) || n != cov.GetLength(1))
                throw new System.IndexOutOfRangeException("Weights and covariance matrix have different dimensions");
            double variance = 0;
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < n; ++j)
                    variance += x[i] * cov[i, j] * x[j];
            
            return Math.Sqrt(variance);
        }

        /// <summary>Calculate expected log return E(ln(R_p))</summary>
        /// <param name="x">Portfolio instrument weights</param>
        /// <param name="r">Instrument expected returns</param>
        /// <param name="cov">Covariance matrix</param>
        /// <returns>POrtfolio expected log return: E(ln(R)) = sum_i x_i * r_i - 1/2 sum_ij x_i * cov_ij * x_j</returns>
        public double calcExpLogReturn(double[] x, double[] r, double[,] cov)
        {
            double ret = this.calcExpReturn(x, r);
            double vol = this.calcVolatility(x, cov);
            return ret - 0.5 * vol * vol;
        }
    }
}
